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Tantalum-oxide-based bi-layered resistance- change memories (RRAMs) have recently improved greatly 
with regard to their memory performances. The formation and rupture of conductive filaments is generally 
known to be the mechanism that underlies resistive switching. The nature of the filament has been studied 
intensively and several phenomenological models have consistently predicted the resistance- change 
behavior. However, a physics-based model that describes a complete bi-layered RRAM structure has not yet 
been demonstrated. Here, a complete electro-thermal resistive switching model based on the finite element 
method is proposed. The migration of oxygen vacancies is simulated by the local temperature and electric 
field derived from carrier continuity and heat equations fully coupled in a 3-D geometry, which considers a 
complete bi-layered structure that includes the top and bottom electrodes. The proposed model accurately 
accounts for the set/reset characteristics, which provides an in-depth understanding of the nature of resistive 
switching. 



Oxide-based resistance-change memories (RRAMs) have been studied intensively during the past decade 
and have attracted increasing interest as potential next-generation nonvolatile memories 1 . To explain 
their resistive switching behaviors, researchers have proposed various mechanisms, among which the 
concept of formation/rupture of conductive filaments (CFs) has been widely accepted 2 . The transition from a 
high- resistance state (HRS) to a low- resistance state (LRS) is called the set process, which is interpreted as a defect- 
induced soft breakdown associated with the migration of non-lattice oxygen ions toward the electrode 3 . In 
contrast, the transition from the LRS to the HRS is called the reset process, in which electric-field-enhanced 
oxygen ion-hopping ruptures the CFs 4 . In the case of bipolar resistive switching, the voltages for the set and reset 
processes have the opposite polarities, thus revealing the role of voltage-driven migration of oxygen ions 5 " 8 . 
Moreover, recent studies have pointed out the crucial role of temperature in the voltage -driven migration process 
at the moment of set or reset transition 910 . 

On the basis of comprehensive and intensive studies, the nature of the filament during the set and reset 
processes is now comprehensible. The presence of one or more localized CFs has also been confirmed by direct 
microscopic methods 1112 and indirect electrical measurements. In addition, several phenomenological or physical 
models have predicted the consistent resistance-change behavior (i.e., current-voltage (I-V) characteristics) of 
RRAMs 5,8 ' 9 . However, existing models have apparently overlooked an important fact. As mentioned earlier, the 
thermal effect is an essential factor of oxygen ion migration during the set and reset processes. Therefore, to 
predict the temperature inside the CF accurately, the thermal properties of both the electrodes and the oxide layer 
should be included in a model. However, most previous models have ignored the effect of electrodes and assumed 
either the CF/electrode interface to be at room temperature 5,9 or the temperature inside the CF to be constant 13,14 . 
Although Menzel et al. have recently proposed an advanced electro -thermal model that considers both the top 
and the bottom electrodes to calculate the temperature of the CF precisely 10 , a physical model that predicts the I-V 
characteristics of an RRAM accurately in a complete bi-layered RRAM structure has not yet been demonstrated. 

In this work, a numerical model, which is based on temperature- and field-accelerated migration of oxygen 
vacancies, is presented for bipolar RRAMs. A 2-D axis -symmetric finite element simulation model, which allows a 
quantitative analysis of field and temperature contributions, accounts accurately for the set/reset characteristics of 
both the DC and the AC responses. The proposed model also provides a microscopic interpretation of inter- 
mediate set/reset states based on differences in the CF morphology. 
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Figure 1 | Schematic diagram and simulated geometry, (a) Schematic 
and cross-sectional TEM images of the Pt/Ta 2 0 5 /TaO x /W bi-layered 
RRAM device, (b) Typical I-V characteristic of Pt/Ta 2 0 5 /TaO x /W device, 
(c) Simulated geometry used in calculation. The axisymmetric geometry 
reduces the problem from one of three dimensions to one of two 
dimensions. 

Results 

The tantalum-oxide-based bi-layered RRAM consisted of a high 
resistive Ta 2 0 5 layer on top of a less resistive TaO x layer (see 
Methods for details), which was sandwiched by a top electrode 
(TE) of Pt and a bottom electrode (BE) of W, as shown in 
Fig. 1(a). When a large negative voltage is applied to the pristine 
device (a process called "forming" shown in Fig. 1(b)), positively 
charged oxygen vacancies are attracted from the TaO x layer to the 
TE, forming a CF 15 ' 16 . This CF is described as a doped region where 
oxygen vacancies (V Q ) act as dopants and thus contribute to local 
electrical and thermal conductivity values 17 . The simulation, which is 
started immediately after the forming process, features a continuous 
CF that connects the top and bottom electrodes (Fig. 1(c)). After the 
forming process, the set and reset processes are described through V G 
migration induced by the local electric field and temperature due to 
Joule heating 517 . Thus, simulation of the set/reset transitions requires 
the self-consistent solution of three partial differential equations 
(PDEs): (1) a drift/diffusion continuity equation for V G transport, 
(2) a current continuity equation for electrical conduction, and (3) a 
Fourier equation for Joule heating, as shown in Fig. 2. 



• Parameters from measurements and assumptions - Eq.(5) 
<7 = cr 0 exp(-E AC I kT) Conductivity [Q- 1 cnr 1 ] 

k lh Thermal conductivity [Wrrr 1 K- 1 ] 





Figure 3 | Parameters from measurements and assumptions. 

(a) Electrical conductivity pre- exponential factor g 0 , (b) assumed 
activation energy for conduction E ACy (c) measured E AC at both HRS and 
LRS, and (c) assumed thermal conductivity k th as a function of local V G 
density n D . 

To describe the drift/diffusion of V G migration, a simple ID rigid 
point ion model by Mott and Gurney is employed 18 . Although the 
ion-hopping model by Mott and Gurney is applicable oxygen ions 
rather than oxygen vacancies, here, we believe that both oxygen 
vacancies and ions can be described by the same equations except 
for the detailed physical parameters. Diffusion is given by D = 1/2 a 2 f 
exp(-EJkT), and the drift velocity is given by v = afexp(-EJkT) 
sinh(qaE/kT), where/is the attempt-to-escape frequency (10 13 Hz) 18 , 
a is the effective hopping distance (0.05 nm), and E a is the activation 
energy of migration. In fact, the effective hopping distance (a) has a 
range of a few nanometers (~1 nm) 18 . However, the value of 
0.05 nm was used in this work for the purpose of best fitting the 
experimental data. This inconsistency is presumably due to differ- 
ences in the physical nature of oxygen vacancies and oxygen ions. 
Considering the drift and diffusion, the time-dependent evolution of 
V G concentration (n D ) can be expressed by the continuity equation, 
as follows: 



• Dependent variables 

n D Concentration of V 0 [cm 3 ] 
T Temperature [K] 
\// Potential [V] 



• Constants 

a Hopping distance, 0.05 nm /? Mesh size, 0.5 nm 

f Attempt-to-escape frequency, 10 13 Hz 

E a Diffusion barrier, 1 .5 eV (ip>0) or 2.5 eV (ujj<0) 

A Coefficient of generation rate, 10 24 cnr 3 s~ 1 



• Oxygen vacancy transport 

Eq.(1) f^L = V-(DV« D -Vrt D ) + G 



• Current continuity 

Eq.(2) V.dV(// = 0 

• Heat (Joule heating) 

Eq.(3) -V-k lh VT = J-E = cr\Vl//\ 2 



Eq.(1) 
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• Parameters from physics- Eqs.(4) 
D = \/2-a 2 -f-exp(-E a /kT) 
v-a- f • Qxp(-E a I kT) • sm\\{qaE I kT) 
G = A Qxp(-(E h - qpE) I kT) 



Diffusivity of V 0 [cm 2 s ~ 1 ] 
Drift velocity of V 0 [cm/s] 
Generation rate of V 0 at iy<0 [cm 3 s " 1 ] 



Figure 2 | The equations and parameters in the proposed model. Three 
PDEs are self- consistently solved within a numerical solver. 
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Figure 4 | Simulated boundary conditions and calculated results. 

(a) Cross section of the simulated cell. A uniform doping concentration of 
n D = IX 10 22 cm" 3 was defined within the CF and TaO x layer as the initial 
forming state, (b) Measured and calculated I-V characteristics. 
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Figure 5 | Simulation results for the reset transition, (a) Calculated map of V Q density (n D ). Calculated profiles of (b) n D , (c) T, and (d) \j/ for states A 
(1.0 V),B(1.7 V),andC(3.0 V) as depicted in Fig. 4(b). The position of z = 15 nm indicates the Ta 2 0 5 /TaO x interface. The depleted gap, defined for n D 
< 5 X 10 21 cm" 3 , is marked by the shaded area. 



_^ =V -(DVn D -vn D ) + G. (1) 

Here, G is the generation term to describe the set process. The kin- 
etics of the set process can be interpreted as a dielectric soft break- 
down associated with the migration or generation of V G through 
thermally activated hopping. The CF growth rate is assumed to 
be proportional to the ion migration rate, which is given by G = 
A exp(-(E b - qVj/kT) = A exp(-(E b - q$E)lkT) l%2 \ where A is a 
pre-exponential constant, E b is the energy barrier for ion hopping 
(1 eV), and /? is the mesh size in the simulation (0.5 nm). Term qfiE 
describes the lowering of the barrier attributed to the applied field. 
This lowering forces the V G to move along the direction of the electric 
field. To express the set process that occurs only at a negative bias, the 
generation term is also applied exclusively at a negative bias. 
Equation (1) can be solved with the current continuity equation 
for electrical conduction, as follows: 

V-(jViA = 0 (2) 

where a is the electrical conductivity, and with the steady- state 
Fourier equation for Joule heating 

-V-^Vr=/-£ = <7|V<A| 2 (3) 

where k th is the thermal conductivity. The transient contribution to 
the Fourier equation was disregarded, given the relatively fast ther- 
mal response of approximately 1 ps (see Supplementary Information) 5,9 . 

The three PDEs in equations (l)-(3) are self- consistently solved by 
a numerical solver (COMSOL) to calculate n D , \jj (potential), and 
T (temperature). To solve equation (3), models for electrical 



conductivity (a) and thermal conductivity (k th ) in the oxide and in 
the CF are required. To this end, both the electrical and the thermal 
conductivities are assumed to depend on n D , as shown in Fig. 3 5 ' 9 . The 
electrical conductivity is given by the Arrhenius equation 21 , a = 
<j 0 exp(— E AC /kT), where <j 0 is a pre-exponential factor and E AC is 
the activation energy for conduction. As shown in Figure 3(a), <j 0 is 
assumed to linearly increase from 5 to 450 Q _1 cm _1 with an increase 
in n D ; these values correspond to the measured values of resistivity at 
both the LRS and the HRS, respectively. In addition, Figure 3(b) 
shows the conduction activation energy E AC used in the calculations. 
The activation energy is 0.016 eV for high n D and linearly increases 
to 0.052 eV with a decrease in n D , consistent with the measured data 
shown in Figure 3(c). Moreover, a linear dependence of k th on n D is 
assumed on the basis of the Wiedemann-Franz law, as shown in 
Figure 3(d) 5 . The minimum value for n D = 0 refers to the thermal 
conductivity of the insulating Ta 2 0 5 , k Ta2 o5 = 0.12 Wm^K" 1 for T 
= 300 K 22 . The maximum k th value at high n D corresponds to that of 
the metallic CF, i.e., thermal conductivity of tantalum k Ta = 
57.5 Wm _1 K _1 . Here, a maximum V G density of 1 X 10 22 cm" 3 is 
chosen approximately where it corresponds to a relative atomic con- 
centration of 20%. Although the particular choice of the maximum 
V G value and the linear approximations of cj 0 , E AC , and k th appear to 
have a weak physical background, the calculation results based on 
these assumptions show good consistency with the experimental 
data; therefore, these assumptions do not limit the validity of the 
calculations. 

Figure 4(a) shows the simulated region of the tantalum-oxide- 
based bi-layered RRAM. In the actual calculations, the axisymmetric 
geometry allowed the 3-D problem to be reduced to a 2-D solution 
with a radial coordinate and a vertical coordinate. The active oxide 
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Figure 6 | Simulation results for the set transition, (a) Calculated map of V G density (n D ). Calculated profiles of (b) n D , (c) E, and (d) G for states 
D (-0.25 V), E (-1.5 V), and F (-0.25 V), as depicted in Fig. 4(b). 



materials (Ta 2 0 5 and TaO x layers) are contacted by two electrodes, 
and all layers including the two electrodes are considered in the 
calculations (Pt TE, o Pt = 10 5 Q^cm" 1 , k Pt = 71.6 Wm^K" 1 ; W 
BE, g w = 2 X 10 5 Q^crrT 1 , k w = 173 Wm^K" 1 ). The boundary 
conditions for equation (2) are \jj = 0 and ijj = V at the bottom and 
top electrodes, respectively. The outermost boundaries of the two 
electrodes are defined as ideal heat sinks with boundary conditions 
T = 300 K; this simplifying assumption is reasonable because the 
electrode area is generally large with respect to the CF. For V G drift/ 
diffusion, no V G flux was assumed at the TE/Ta 2 0 5 and TaO x /BE 
interfaces. In addition, only diffusion was considered in the TaO x 
layer (i.e., the drift term (v) of equation (1) is zero) because this layer 
only serves the role of an oxygen vacancy reservoir 1516 . A uniform 
doping concentration of n D = 1 X 10 22 cm" 3 was defined within the 
CF and TaO x layer as the initial state (i.e., the forming state). The CF 
size was set to a diameter of 10 nm, in agreement with direct evalua- 
tions by using conductive atomic force microscopy 23 . 

Discussion 

Figure 4(b) shows the measured and calculated I- V characteristics 
during the set and reset processes after the forming process. The reset 
transition starts at 1.85 V and the resistance gradually increases, 
finally achieving a difference of roughly one decade after V RESET = 
3.0 V. Similarly, the set transition starting at — 0.5 V also shows good 
consistency with the measured data. Figure 5 shows calculated (a) 
2-D maps of n D as well as 1-D profiles of (b) n D , (c) T, and (d) \jj at 
bias points A, B, and C. As the voltage increases during the reset 
sweep, the temperature increases as a consequence of Joule heating, 
as shown in Fig. 5(c). V G migration is localized at the point of max- 
imum temperature (i.e., near the TE) owing to the activation of drift 



velocity by exponential T. Consequently, V G migrates toward the BE, 
resulting in the opening of a depleted gap. The opposite process 
occurs in the set transition. Figure 6 shows the calculated (a) 2-D 
maps of n D and the 1-D profiles of (b) n D , (c) E (electric field), and (d) 
G (generation rate of V G ) at the bias points D, E, and F. Given that the 
electric field at the depleted gap increases with an increase in the 
negative bias, the generation rate of V G increases, as shown in 
Figs. 6(c) and 6(d). Consequently, the depleted gap at the CF is re- 
filled by V G , after which the high electric field is no longer applied to 
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Figure 7 | Temperature dependency, (a) Calculated I-V characteristics at 
different device temperatures. Measured and calculated data for the (b) 
reset transition voltage and (c) resistance value of HRS. 
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Figure 8 | Gradual reset process, (a) Measured and (b) calculated I-V 
curves, showing reset transition, obtained by varying V RESET (1.8 V to 
3 V). 

the CF. This suppresses further V G generation or migration and is 
responsible for a self-limiting set transition. 

By applying different boundary conditions of temperature, the 
thermal effect on the I- V characteristics can be predicted by calcula- 
tions. Figure 7(a) shows the calculated I-V curves under different 
temperature conditions. The reset transition voltage (i.e., the 
moment at which the resistance starts to decrease) noticeably 
decreases with an increase in the device temperature. In addition, 
the resistance at the HRS increases with an increase in the device 
temperature. The reset transition occurs at a lower bias voltage 
because V G migration (drift and diffusion) is more favorable at higher 
device temperatures. Moreover, the difference of conduction curve at 
the HRS is due to the temperature-dependent conductivity enhance- 
ment of a through the Arrhenius equation. These two distinctive 
features agree well with the experimental data shown in Figs. 7(b) 
and 7(c), where the measurements are carried out in the range of 
— 20°C-85°C; accordingly, the simulation accurately reflects the 
thermal effect for the set and reset transitions. To further verify the 
simulated results, Figures 8(a) and 8(b) show the measured and 
calculated I-V curves, including reset sweeps stopped at various 
values of Vreset. The simulation accurately reproduces the different 
features of resistive switching with different values of V RESET and 
predicts a gradual reset transition behavior. 
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Figure 10 | Calculated I-V curves for filaments of different shapes. 

(a) rectangular, (b) trapezoidal, and (c) inverted trapezoidal. Despite the 
different shapes, the position of the depleted gap is almost identical. In 
addition, the I-V" characteristics are also independent of the filament shape. 

The application of the proposed model is not limited to repro- 
ducing the DC I- V characteristics. Time-dependent reset character- 
istics were further explored by investigating the response to square 
pulses. Time-resolved experiments were carried out in which a series 
of set/reset pulses was applied to the RRAM device with a serially 
connected resistor, and the response to an AC pulse was simulta- 
neously detected by an oscilloscope. Figures 9(a) and 9(b) show the 
calculated and measured response profiles to the applied reset pulse 
as a function of time. When a single square pulse for the reset trans- 
ition is applied to the device, the response current increases following 
a voltage increment because the initial state of the device is in the 
LRS. Over time, the response current rapidly decreases, which indi- 
cates the moment of reset transition. The measured time to the reset 
transition (t RES ET) is 1.3 jas for the input pulse shown in Fig. 9(b), 
which is comparable to the calculated value of t reset = 3.5 |is. 
Consequently, the simulation can describe the AC resistance-change 
response reliably in addition to the DC I-V characteristics. 

Although the morphology of the CF has been observed by direct 
imaging with transmission electron microscopy (TEM), the effort 
required was extensive and the extracted results were insufficiently 
informative for comprehensively understanding the mechanism of 
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Figure 9 | Transient response, (a) Measured and (b) calculated AC 
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Figure 11 | Effect of different TE materials, (a) Calculated I-V curves for 
different TE materials, (b) 2-D maps of n D corresponding to TE materials, 
and temperature profiles in the CF and TaO x layer. 
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resistive switching 11 . In contrast, a merit of the simulation is that it 
provides a microscopic interpretation of intermediate set/reset states 
based on differences in the CF morphology. As mentioned above, the 
thermal properties of both the electrodes and the oxide layer should 
be included in the model to accurately describe the resistive switch- 
ing phenomenon. Previously reported models predicted that the 
depleted gap at the reset transition should be positioned in the mid- 
dle of the CF when temperature- driven process is dominant 4,917 . On 
the other hand, several papers have verified by electrical measure- 
ments that the depleted gap is located near the electrode/oxide inter- 
face. To overcome this disagreement between the model and the 
measurement result, modified models have been explained that the 
depleted gap can be located near the electrode/oxide interface, par- 
ticularly when the shape of the CF is conical 24 ' 25 . However, the ana- 
lyses based on these modified pre-existing models are also inaccurate 
because they did not consider the thermal property of the electrodes. 
Figure 10 shows the I- V characteristics and 2-D maps of n D corres- 
ponding to filaments differing in their initial shape: (a) rectangular, 
(b) trapezoidal, or (c) inverted trapezoidal. It should be noted that the 
location of the depleted gap in the CF is independent of the shape. In 
particular, the filament with an inverted trapezoidal shape changes to 
a rectangular one after consecutive switching cycles. These results 
imply either that the depleted gap is not always located in the middle 
of the CF or that its location depends on the shape of the filament. 
Indeed, a crucial factor for the reset transition is the thermal con- 
ductivity ratio between the TE and the oxide. Figure 1 1(a) shows the 
calculated I-V characteristics based on different TE materials, and 
Figure 11(b) shows the corresponding 2-D maps of n D and the tem- 
perature profiles. The position of the depleted gap in the CF is deter- 
mined by the thermal conductivity of the TE, and the DC I-V 
characteristics are also distinct. (The BE material has a negligible 
effect owing to the presence of a relatively thick TaO x layer. In the 
case of a single-layered RRAM, the effect of the BE is also significant). 
As the thermal conductivity of the TE increases, the generated heat 
by Joule heating is confined within the oxide layer; consequently, the 
depletion gap of the CF is positioned away from the TE/oxide inter- 
face. Therefore, the effect of the electrodes should be considered in an 
analytical model to predict the nature of resistive switching accur- 
ately. The model proposed here provides an accurate interpretation 
of intermediate set/reset states based on differences in the CF 
morphology. 

In summary, a physical and numerical model has been proposed 
on the basis of temperature/field- driven migration of oxygen vacan- 
cies in bi-layered bipolar RRAM devices. The proposed model allows 
the calculation of resistance- change characteristics for both the DC 
and the AC input signals. By considering the top and bottom elec- 
trodes, a more accurate temperature profile than that of the pre- 
existing models is obtained, which in turn leads to more accurate 
simulation results. The proposed model is potentially useful for gain- 
ing physical insights into the morphology of the CF during set and 
reset transitions, which is valuable for exploring the scaling and 
multilevel capabilities of RRAM devices. 

Methods 

All measured samples (Pt/Ta 2 0 5 /TaO x /W devices) had a 1 um X 1 urn crossbar 
structure. A 50-nm-thick W layer was deposited by sputtering and then patterned by 
dry etching to form the bottom electrode (BE) on a 3-nm-thick adhesion Ti/500-nm- 
thick Si0 2 /Si substrate. Then, a 20-nm-thick conducting TaOx layer was deposited on 
the BE by DC magnetron sputtering in a reactive 0 2 + Ar ambient, with Ta as the 
target. Subsequently, an insulating Ta 2 0 5 layer was formed by a pulsed 0 2 plasma 
oxidation process that involved a 2-s plasma-on step and a 4-s plasma-off (purge) step 
in a single cycle. The total oxidation time was controlled by the number of cycles, 
which was 80 in the present study. Finally, a 20-nm-thick Pt layer was evaporated and 
patterned by a lift-off process. 

The DC electrical characterizations were carried out by an Agilent B1500A semi- 
conductor analyzer in the voltage- sweep mode. During AC pulse measurements, a 
series resistor ( 1 kQ) was serially connected to the RRAM device. Input pulses were 
generated by a B1500A, and the response was detected by an oscilloscope (Tektronix 
DPO 70604). 
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